Billiards in magnetic fields: A molecular dynamics approach 
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We present a computational scheme based on classical molecular dynamics to study chaotic bil- 
liards in static external magnetic fields. The method allows to treat arbitrary geometries and 
several interacting particles. We test the scheme for rectangular single-particle billiards in mag- 
netic fields and find a sequence of regularity islands at integer aspect ratios. In the case of two 
Coulomb-interacting particles the dynamics is dominated by chaotic behavior. However, signatures 
of quasiperiodicity can be identified at weak interactions, as well as regular trajectories at strong 
magnetic fields. Our scheme provides a promising tool to monitor the classical limit of many-electron 
semiconductor nanostructures and transport systems up to high magnetic fields. 

PACS numbers: 05.45.Pq,82.40.Bj,73.21.La 
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I. INTRODUCTION 

Classical and quantum billiard systems [H, 0] are of 
significant interest both in nonlinear physics and in ap- 
plications based on low-dimensional nanostructures [3(. 
For example, quasi-two-dimensional (quasi-2D) quantum 
dots [H are studied in view of emerging applications in 
the field of quantum computation [5j]. They exhibit de- 
terministic ballistic motion of the electrons as "billiard 
balls" and provide the possibility to tune their shape, 
size, and electron number. A particularly intriguing fea- 
ture is the connection between classical dynamics and 
the statistical properties of the corresponding quantum 
system [|| 0] ■ For systems with mixed chaotic and reg- 
ular dynamics, the Berry-Robnik formula Q links the 
volume ratio of regular and chaotic regions in classical 
phase space to the quantum-mechanical level distribu- 
tion [I El]. 

External magnetic fields pose, on the one hand, an 
interes ting complication to classical (and quantum) bil- 
liards [111 ], and, on the other hand, provide an eas- 
ily accessible way to experimentally control the par- 
ticle dynamics. Recently, magnetic fields have been 
used to manipulate electron transport in coupled elec- 
tron billiards LUJ . In many cases^ 
or triangular [lj 



lar \m UJ, U5 



e.g., in rectangu- 
billiards, an ex- 



ternal magnetic field leads to mixed dynamics between 
regularity and chaoticity. The breaking of time reversal 
symmetry due to the presence of a magnetic field results 
in new properties of the level spa cing statistics of the 
corresponding quantum system [18|, Qj| ■ 

In contrast to freely tunable parameters, such as ex- 
ternal magnetic and electric fields, interactions between 
particles are inevitably present in any realistic physical 
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system. While single-particle billiards have been stud- 
ied thoroughly for many years now, billiards of interact- 
ing particles are still a relatively young field. Classical 
billiards for two interacting particles have been studied 
using various models, ej^ Coulomb-like interactions in 
a one-dimensional box (20| and in an isotropic [2l| and 
anisotropic harmonic oscillator [22j, as well as apply- 
ing hard-sphere contact interaction in a rectangular [231 ] 
and a mushroom-shaped box (24|. The statistical me- 
chanics of such systems has also been extensively stud- 
ied recently [25| . Quantum-mechanically, interaction- 
induced chaos has been studied in a two-electron quan- 
tum dot [21], [22], HH , and, very recently, also in the frame- 
work of time-dependent density- functional theory [27l.[28j 
- an approach that might enable examination of quantum 
chaos in systems containing a large number of interacting 
particles. 

Single-particle billiards have traditionally been studied 
by either reducing the dynamics of the system to a bounc- 
ing map (for magnetic single-particle billiards, see, e.g., 
[29|), or by investigating the infinitesimal variations of 
the trajectories using the method of Jacobi fields [30l. [3~TI] . 
In an interacting billiard, however, the trajectory of a 
particle between successive bounces is not known in ad- 
vance, as its motion is coupled to the motion of all other 
particles. The locations of the bounces at the wall are not 
given by simple geometric considerations anymore, and 
thus the methods used to study single-particle systems 
do not carry over in a straightforward way. 

In this paper, we present a classical molecular dynam- 
ics scheme that allows to calculate the trajectories of in- 
teracting particles in an arbitrary 2D billiard system ex- 
posed to a uniform and perpendicular magnetic field. To 
demonstrate the method, we focus on single- and two- 
particle dynamics in rectangular billiards. In the single- 
particle case, we present an efficient method to system- 
atically obtain "regular" and "chaotic" regions in phase 
space, which allows us to monitor the combined effect 
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of the magnetic field and the rectangle shape. We find 
a pattern of increased regularity at integer aspect ratios. 
In the two-particle case mostly chaotic behavior is found, 
but also regular orbits at high magnetic fields. The rel- 
evance of the method in studying the classical limit of 
collective effects in many-electron structures is discussed. 

II. METHOD 
A. Propagation of particles 

To calculate the trajectories of charged particles, we 
use a modified velocity vcrlct algorithm suited for incor- 
porating arbitrarily strong static homogeneous external 
magnetic fields [33|. With a magnetic field B = (0, 0, B) 



pointing in z direction, the acceleration of a charged par- 
ticle reads 

a(t) = a c (()-(le 2 xv(f), (1) 

where a if) is the velocity-independent part of the ac- 
celeration depending only on external forces, and Q = 
qB/m is the cyclotron frequency for a particle with 
charge q and mass m. We use Hartree atomic units 
throughout the paper, such that H = e = m e = 
l/(47reo) = 1 and the velocity of light has the value 
c as 137.036. Furthermore, the factor 1/c in the Lorentz 
force law is absorbed into B, such that we have Q, = B 
for electrons. Within the modified velocity verlet algo- 
rithm presented in Ref. [3^ . each particle is propagated 
using the following equations: 
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B. Phase space maps for single-particle billiards 

In a single-particle billiard system, the kinetic energy, 
and consequently also the velocity v = iv^+Vy) 1 / 2 of the 
particle, is a constant of motion. The dynamics of the 



billiard is determined by the boundary conditions (see 
below) and the relative strength of the magnetic field. 
The latter quantity is here given by a parameter 



fx — R c /L x , 



(7) 
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where R c = v/B is the cyclotron radius and L x is 
the length of one side of the system (here a rectan- 
gle). The constant of motion can be used to reduce 
the four-dimensional phase space (x, y, v x , v y ) to a three- 
dimensional (3D) one, where we have chosen the space 
spanned by (x, y, v x )- To identify regular and chaotic re- 
gions in this phase space, we use the following procedure: 

1. We choose a 2D cross section (x,v x ) through the 3D 
phase space and divide it into a number of cells. 

2. For one cell, we pick two phase space points in the 
cell that are very close to each other, but not iden- 
tical up to the numerical precision. 

3. We follow the trajectories through these two 
points for a certain propagation distance stot us- 
ing Eqs. ([I])-©, and record all cells in the cross 
section through which they pass. 

4. After having propagated for a distance s t ot, we cal- 
culate the distance between the points in phase 
space, which is a measure of the "regularity" of 
the trajectory. We save this distance to all cells 
we have passed. If a cell has already been passed 
by a previous run, we take the maximum of the 
distances. 

5. We start over from point (2) by picking another cell 
that has not yet been traversed by a trajectory, and 
repeat the whole process until all cells have been hit 
by a trajectory at least once. 

6. We then plot the distances stored in the cells of 
our 2D cross section as a color-coded "matrix plot" . 
In the following, these plots will be called "phase- 
space maps" . Small numbers correspond to "regu- 
lar" phase space cells, large numbers to "chaotic" 
cells (see below for details). 

The algorithm can be efficiently parallelized, because 
trajectories originating from different cells can be prop- 
agated independent of each other. Our code uses the 
Message Passing Interface (MPI) and a master/slave 
paradigm. The master process keeps track of the phase 
space map and distributes free cells, i.e., cells that have 
not yet been hit by any trajectory) to the workers. The 
workers perform the propagation of the trajectories and 
communicate the traversed cells and the phase space dis- 
tance after the propagation distance Stot back to the 
server. 



III. RESULTS 

A. Single particle 

We demonstrate our computational scheme by con- 
sidering rectangular billiards with side lengths L x = 1 
(fixed) and L y = f3 L x (varied) , where (3 is the aspect ra- 
tio. The strength of the external magnetic field has been 




p (%) 



FIG. 1: (color online), (a) Phase space map for a rectangular 
billiard with aspect ratio (3 = 2. The color map indicates 
the phase-space distance A between two orbits having a small 
initial perturbation, (b) Ordered phase-space distances for all 
the cells plotted in (a). The dashed line shows the threshold 
(A = 0.1) between chaotic and regular motion. 



fixed to B = 1, so that the cyclotron radius R c = v/B 
is determined by varying the velocity of the particle. In 
the single-particle case, we focus on the dynamics of the 
system as a function of (3 and [i = R c /L x . In both of 
the limits fx — > and fx — > oo the motion is regular, 
the former corresponding to infinitely many circular or- 
bits (cf. Landau-level condensation in confined quantum 
systems) and the latter corresponding to linear motion at 
zero field, which is always regular in rectangular billiards. 
At < /i < oo the dynamics is generally mixed except 
at particular values of jj, when the system is completely 
chaotic 

Fig. []Ja) shows an example of a phase space map cal- 
culated for the parameters (3 = 2 and \x = 0.6. The 
scheme described in Sec. Ill B I has been used to calculate 
the figure. The cross section through the phase space 
has been partitioned into 150 cells in each direction (x 
and v x )- The color scale indicates the phase-space dis- 
tance A after propagating the trajectories by a distance 
of stot = 60 

We find distinct areas of regularity associated with 
KAM (Kolmogorov-Arnold-Moscr) islands [l|. Fig- 
ure [IJb) shows the phase-space distances of all cells 
sorted in ascending order. The sharp onset of the 
curve indicates a distinct separation between regular and 
chaotic motion. To consistently determine this separa- 
tion, we choose a threshold of A = 0.1 shown in the 
figure as a dashed line. Thereby, this particular system 
is regular by a fraction of 25%. We presume that by using 
a very high resolution it should be possible to determine 
and c ateg orize phase-space cells corresponding to weak 
chaos [33] ]. This topic is, however, beyond the scope of 
this work and left for future research. 

In Fig. [5] we show the proportions of regularity, es- 
timated as shown in the example in Fig. [TJ for square 
(/3 = 1) billiards as a function of /j,. We find excellent 
agreement with the result of Berglund and Kunz [l3| that 
has been calculated using an exact method. This con- 
firms the accuracy of the proposed scheme up to strongly 
curvilinear motion, i.e., small values of fi. Hence, we ex- 
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FIG. 2: Proportions of regularity in square billiards as a func- 
tion of /i = R c /L x , i.e., the ratio between the cyclotron radius 
and the side length. 
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FIG. 4: (color online). Initial configuration of the calcula- 
tions for two Coulomb- interacting particles in square billiards 
subjected to a perpendicular magnetic field. 



B. Two particles 
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FIG. 3: (color online) Proportions of regularity in rectangular 
billiards as a function of fi and the aspect ratio (3. 



pect the method to be reliable also in more complicated 
systems with many particles and/or different boundaries. 

To assess the effect of the billiard shape onto the dy- 
namics, we have calculated the proportions of regularity 
as a function of both /i and the aspect ratio (3. The re- 
sult is shown in Fig. [3] for 0.3 < (3 < 3.2 in steps of 0.1 
and for 0.2 < \i < 2 in steps of 0.01. Fig. [3] required the 
calculation of 5430 phase space maps, each consisting of 
22 500 cells, thereby demonstrating the numerical effi- 
ciency of the scheme. Note that Fig. [3] is not symmetric 
around (3 = 1, because we have varied L y = (3 and thus 
the system area is not kept constant. 

Wc find several islands of increased regularity centered 
at f3 = n/2 with n = 1,2,3,.... Overall, the "most 
regular" case is the square billiard ((3 = 1), as expected. 
A more detailed analysis of the regularity patterns and 
their connections to the periodic orbits will be performed 
elsewhere. 



We now turn to the dynamics of two particles inter- 
acting via Coulomb repulsion in a square well (/3 = 1). 
Now the velocities (and thus also the cyclotron radii) 
are no longer constants of motion. The phase space is 
eight-dimensional, and instead of the phase space map 
described in section Hi B[ we calculate so-called bouncing 
maps by recording the values (x, v x ) corresponding to the 
bounces of one of the particles on the lower boundary 
(y = 0) of the system. 

We investigate the dynamics with different values for 
the ratio 
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fik(* = 0) 
E p (t = 0) 



(8) 



for the initial configuration, where E p = |r — r'| _1 /2 is 
the Coulomb potential energy and = v 2 /2 is the ki- 
netic energy. The quantity 7 essentially determines how 
"strongly interacting" the system is, as it fixes the aver- 
age ratio of Ek(t) and E p (t) for the full time- dependent 
system through the initial energy components. In phys- 
ical applications this ratio could be varied by changing 
either the particle density or the system size. A well- 
known example of the limit where the potential energy 
dominates is the Wigner crystal [34[ forming in the elec- 
tron gas at low densities. 

In the following examples we have fixed the initial 
positions of the particles to {x\,yi) = (2/5,3/10) and 
(^2,2/2) = (7/10,7/10). The initial Coulomb energy in 
this case is E p (t = 0) = 2. After fixing 7 in Eq. ©, 
the initial kinetic energy Ek{t = 0) is distributed equally 
to both particles, and the initial velocities Vi = V2 = 
(0, \/Ek) point in the y direction. The initial configura- 
tion is visualized in Fig. [U 

The remaining parameter to be fixed is fi(t = 0) de- 
fined in Eq. ([7]). Note that again we fix only the initial 
condition, and in the time-dependent run, the values of 
fi for both particles vary due to changes in the velocities. 
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FIG. 5: (color online). Upper panel: Classical trajectories for 
two relatively weakly interacting particles indicated by blue 
(black) and red (gray) colors in square billiards. The magnetic 
field is zero. Lower panel: Bouncing map for the particle with 
the blue (black) trajectory in the upper panel. 

Since the initial velocity is determined through 7, we fix 
fi(t = 0) through B in contrast with the single-particle 
case where we always had B = 1. 

First, we set 7 = 30 and the magnetic field to zero 
(/i — > 00) and propagate sufficiently long to obtain a 
bouncing map with a large number of points. Figure [5] 
shows the trajectories of the particles up to t = 5 (up- 
per panel) and the bouncing map up to t = 3 x 10 4 
(lower panel). The number of bounces is ~ 7.8 x 10 5 . 
Apart from a few exceptions, the particles remain sep- 
arated in the left and right parts of the system due to 
the Coulomb repulsion. However, as the interaction is 
relatively weak, both particles move in the y direction, 
almost undisturbed from their initial conditions. Close 
to the left and right boundaries, where the interaction 
is weakest, the dynamics is most regular. This can be 
seen in the trajectories, which are almost straight lines 
in that regime. In addition, the bouncing map shows reg- 
ular curvilinear albeit blurry zones (see the inset in the 
lower panel of Fig. [3]). These features may be designated 
as quasi- regular motion in the system [33j. 



FIG. 6: (color online). Same as Fig. [5] but for relatively 
strongly interacting particles. 



In the following example we keep the magnetic field at 
zero but increase the relative amount of interaction en- 
ergy such that 7 = 1/30. The trajectories and bouncing 
map are shown in Fig. [6l Here, the dynamics is very dif- 
ferent from the weakly interacting case. Both particles 
occupy the whole area of the system, but due to their 
strong repulsion, the corners are considerably more oc- 
cupied than the central region, which is characterized by 
"scattering" trajectories of high curvature. The bouncing 
map in the lower panel is completely chaotic. Increasing 
the interaction even further would enable to study clas- 
sical Wigner crystallization (34j in a dynamic picture. 
In the present system, for example, the Wigner crystal 
would consist of two diagonal configurations summed up 
to a four-point crystal. 

Finally, we consider two systems with 7 = 2, where 
the magnetic field is set to values corresponding to ji(t = 
0) = 1/4 and fj,(t = 0) = 1/32, respectively. The tra- 
jectories are plotted in Fig. [7] In the first case (a) the 
system seems to be fully chaotic, whereas the latter con- 
figuration (b) leads to regular isolated orbits forming a 
ring-like structure. In this case, the "interaction axis" 
(i.e., the dashed line in Fig. [?]) performs a circular mo- 
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FIG. 7: (color online). Classical trajectories for two interact- 
ing particles at two different magnetic fields corresponding to 
p = 1/4 and /i = 1/32, respectively. 

tion that is superimposed by strongly confined cyclotron 
motions at the opposite ends of the axis. The character- 
istics of this motion are further illustrated in the inset 
of Fig. [7Jb), which shows the trajectories soon after the 
beginning of the time propagation. The — at a first 
glance — counterintuitive result that a repulsive interac- 
tion leads to bound motion can be understood by consid- 
ering the combined effect of Coulomb repulsion and the 
strong magnetic confinement through the cyclotron mo- 
tion. When the particles increase their relative distance, 



the gain in kinetic energy (at the expense of Coulomb 
energy) results in an increased radius of the cyclotron 
motion. The different curvature of the trajectory on the 
different sides of the "circle" gives rise to a bent cycloidal 
motion which can, for the right choice of the parameters, 
lead to a bound motion as depicted in Fig. [7£b). 

The above results on the classical dynamics in mag- 
netic fields suggest to study the relation to the corre- 
sponding quantum mechanical situation in semiconduc- 
tor quantum dots. In fact, interesting vortex patterns 
and edge localization have been found in rectangular 
many-electron quantum dots at high magnetic fields [H| . 
Such patterns may exist - in a statistical picture - also in 
a classical system. A particularly interesting case would 
be the quantum-mechanical analog of the bound motion 
shown in Fig. [TJb). Moreover, our scheme would allow 
to study the effects of interactions on the classical limit 
of electron transport in billiard arrays fl2| . 



IV. SUMMARY 

We have introduced a computational scheme based on 
molecular dynamics to study classical billiards of inter- 
acting particles in external magnetic fields. The accuracy 
and efficiency of the method has been demonstrated in 
rectangular billiards. We have found excellent agreement 
with numerically exact method in single-particle square 
billiards as a function of the magnetic field. Changing 
the aspect ratio (3 of the rectangle leads to islands of 
increased regularity at j3 = n/2 with n = 1, 2, 3, . . .. In 
square billiards of two interacting particles we have found 
signatures of quasipcriodic orbits at weak interactions, 
and localization at strong interactions. Large magnetic 
fields may lead to regular patterns also for interacting 
particles. The scheme opens up the path to study the 
classical limit of realistic many-particle systems related 
with, e.g., electronic transport experiments in mesoscopic 
structures. 
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